Random walks near Rokhsar-Kivelson points 
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There is a class of quantum Hamiltonians known as Rokhsar-Kivelson(RK)-Hamiltonians for 
which static ground state properties can be obtained by evaluating thermal expectation values for 
classical models. The ground state of an RK-Hamiltonian is known explicitly, and its dynamical 
properties can be obtained by performing a classical Monte Carlo simulation. We discuss the 
details of a Diffusion Monte Carlo method that is a good tool for studying statics and dynamics 
of perturbed RK-Hamiltonians without time discretization errors. As a general result we point 
out that the relation between the quantum dynamics and classical Monte Carlo simulations for 
RK-Hamiltonians follows from the known fact that the imaginary-time evolution operator that 
describes optimal importance sampling, in which the exact ground state is used as guiding function, 
is Markovian. Thus quantum dynamics can be studied by a classical Monte Carlo simulation for 
any Hamiltonian that is free of the sign problem provided its ground state is known explicitly. 



INTRODUCTION 

Microscopic models that describe electronic behavior 
of materials are seldom exactly solvable. Nevertheless it 
is often found that lots of insight can be gained by study- 
ing a solvable model that resembles the accurate micro- 
scopic model in question. A typical approach is to ne- 
glect certain couplings in the accurate microscopic model 
in order to obtain a solvable model. This brings up the 
question of how important the neglected couplings are. 
One way of finding out is to simulate the accurate model 
numerically and compare to the results of the solvable 
model. If the results are sufficiently similar one can, with 
reasonable confidence, claim to understand the physical 
behavior using insights gotten from the solvable model. 

The number of exactly solvable models are however 
relatively small. They tend to be either models of non- 
interacting particles, like free fermions, or non-trivial, but 
one-dimensional. Fortunately understanding and insight 
isn't necessarily tied to exactly solvable models. Models 
that can be fully or partially mapped onto other well- 
known models are also useful as older established results 
can be recycled. 

There is a class of quantum models for which static 
ground state properties can be calculated by evaluating 
thermal expectation values of a corresponding classical 
model having the same number of space dimensions. This 
is a tremendous simplification as much are known about 
classical statistical mechanics models. 

The most well-known of these quantum models is 
the quantum dimer model(QDM) at a special point in 
parameter space; the Rokhsar-Kivelson (RK) point I']. 
There are also many other models with this property 2, 
0, with Hamiltonians known as (generalized) RK- 
Hamiltonians. Beside static properties it is also possible 
to obtain information about excited states and dynami- 
cal properties for RK-Hamiltonians, although in a rather 
unorthodox way. Henley'si showed that dynamical cor- 



relation functions at the RK-point can be gotten by per- 
forming a continuous-time Monte Carlo simulation of the 
classical model using appropriate Monte Carlo dynamics, 
and interpreting time- correlation functions of the clas- 
sical simulation as imaginary-time correlation functions 
of the quantum system. This observation was recently 
utilized in Ref. |6j to determine excitations of the quan- 
tum dimer model on the triangular lattice at the RK- 
point. The observation that quantum dynamics can be 
gotten from classical Monte Carlo simulations is however 
not peculiar to RK-Hamiltonians. It is far more general, 
but requires knowledge of the exact ground state wave 
function. This was recently pointed out in Ref. but 
follows also, as will be explained here, as a consequence 
using importance sampling with an optimal guiding func- 
tion in Quantum Monte Carlo. 

Assuming that the physics at the RK-point is under- 
stood or at least that it can be calculated relatively eas- 
ily, it is interesting to ask if the same physics also holds 
away from the RK-point. An obvious approach to at- 
tempt answering this is to perform a finite-temperature 
path-integral Monte Carlo simulation of the full quan- 
tum model. This is more complicated than perform- 
ing a classical Monte Carlo simulation at the RK-point, 
but can be carried out for most quantum models pro- 
vided the temperature is high enough. This approach 
was followed by Moessner and Sondhi in order to show 
that the spin liquid state existing at the RK-point of the 
QDM on the triangular lattice also extends to other val- 
ues of the parameters However it is difficult to push 
these methods to low temperatures. The non-local up- 
date techniques usually employed to speed up quantum 
Monte Carlo simulations such as loops 8], worms 9] or 
directed-loopsjl^ are not easily applicable to the quan- 
tum dimer model as the configurations on each time-slice 
are heavily constrained. The directed-loop method have 
however been used to study the classical dimer model, 
for which very impressive system sizes can be studied 
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effectively llj. On the other hand the special Monte 
Carlo dynamics implied by the non-local moves doesn't 
correspond to the dynamics of the quantum dimer model 
at the RK-point. 

In this review we will discuss another Monte Carlo 
technique which is an excellent tool to study models 
in the vicinity of RK-Hamiltonians at zero temperature. 
This technique, which is a continuous-time formulation of 
the well known Diffusion Monte Carlo(DMC) technique, 
becomes equivalent to a continuous-time classical Monte 
Carlo method for RK-Hamiltonians. In fact, as is true 
for any DMC method, it becomes equivalent to a classical 
Monte Carlo method whenever it is used in conjunction 
with importance sampling that employs the exact ground 
state wave function as guiding function. 

We concentrate on systems close to RK-points here. 
As is well known with DMC-methods the quality of the 
results depend to a large extent on the quality of the guid- 
ing function. The advantage of being close to a RK-point 
is that the ground state wave function is known exactly 
at the RK-point and can be used as good approximate 
guiding function in its vicinity. 

This review is divided into two main parts. In Sec. „ 
the numerical method is discussed in details. Particular 
emphasis is put on how to extract observables reliably. 
In Sec. B the apphcation of the method to perturbed RK- 
Hamiltonians is discussed, and results pertaining to the 
quantum dimer model on a square lattice are given. 

QUANTUM MONTE CARLO 

There are two main branches of Quantum Monte Carlo 
methods. The first type of methods uses a stochastic pro- 
cess to sample the finite temperature quantum partition 
function and extract observables from this. The main 
challenge using these methods is to engineer efficient up- 
dates as the objects, or variables, in the partition func- 
tion are extended objects, or paths, as follows from the 
path integral formulation of quantum mechanics. A big 
step forward in finding efficient update schemes has been 
achieved using avatars of the Swendsen-Wang0 cluster 
update, so called loopls*!, worm 9] or directed-loop|0| up- 
dates. These Monte Carlo methods are exact in the sense 
that they have no systematic bias of any sort, they are 
easily formulated in continuous-time, and are easily pro- 
grammable. A drawback, but also a source of flexibility, 
is the different ways one can construct these non-local 
updates. This needs to be reconsidered for all models in 
question, although it can be automated to a great ex- 
tent and rules of thumb for how to choose good rules 
exists [HQ. The finite-temperature methods are most 
efficient at high temperatures, and the zero temperature 
behavior is only obtained asymptotically by performing 
simulations at decreasingly lower temperatures. 

The other branch of Monte Carlo methods use stochas- 



ticity to simulate the outcome of (repeated) matrix mul- 
tiplications. This general class of methods is known as 
Projector Monte Carlo(PMC), -repeated multiplications 
of the same matrix projects out the eigenstate with the 
highest eigenvalue. For quantum systems one is inter- 
ested in the lowest eigenvalue of the Hamiltonian. Thus 
in Green function Monte Carlo which is a particular pro- 
jector Monte Carlo technique the iterated matrix is the 
inverse of the Hamiltonian. Another method, DMC, uses 
the matrix exp(— i/r) where H is the Hamiltonian and r 
is imaginary-time. 

For long times, or equivalently many matrix multipli- 
cations, any initial state with some overlap with the true 
ground state will sample the components of the ground 
state wave function and ground state observables can be 
obtained. Moreover because DMC evolves the wave func- 
tion identically to the imaginary-time evolution quantum 
dynamic observables can be obtained easily from the evo- 
lution of the state. 

PMC techniques complement the finite-temperature 
techniques because they are genuine zero-temperature 
methods, -no extrapolations to zero temperature are nec- 
essary. They are also excellent at finding the ground 
state energy in different sectors of conserved quantum 
numbers, something which generally is difficult with 
the methods that work in the grand-canonical ensem- 
ble. Moreover no ingenuity in finding efficient updates is 
needed as everything is directly dictated by the Hamil- 
tonian itself. 

The evolution matrix in DMC is in general not a 
Markovian matrix. Thus DMC cannot be directly inter- 
preted as a classical Monte Carlo method. The standard 
wavjlSj to cope with the lack of probability conservation 
is to include additional branching processes. Unfortu- 
nately these have a tendency to make simulations unsta- 
ble and some feedback control is needed. This feedback 
results in a systematic bias0| to the observables which 
can be removed by reweighting the simulation at the ex- 
pense of additional statistical errors. 

It is often stated repeatedly in the literature that DMC 
cannot be formulated in continuous (imaginary) time, 
and that repeated runs with decreasingly smaller time- 
intervals must be performed in order to quantify the error 
induced by a finite time-step. However the construction 
and implementation of a continuous-time formulation of 
DMC is straight-forward for lattice mo dels [l^. a fact we 
aim to explain in the next section. Other lattice PMC 
techniques^] using different operators to project onto 
the ground state, such as for instance 1 — (i7 — -E'o)''', 
are also free of time-discretization errors. However these 
methods do not have the advantage, as do DMC, that 
quantum dynamics can be obtained easily. For another 
newly developed continuous-time Monte Carlo method, 
see Ref. [l9| . 
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Continuous-time Diffusion Monte Carlo algorithm 

The basic idea behind DMC is to simulate the power 
method stochastically. The power method uses repeated 
matrix multiplications to project out the eigenstate hav- 
ing the largest eigenvalue. In DMC the imaginary-time 
evolution operator projects out the ground state of the 
Hamiltonian. Specifically 



(1) 



where H is the Hamiltonian, is the ground state en- 
ergy and \xi) is an arbitrary initial state having overlap 
Cxi = {'4'a\xi) with the ground state |'(/'o)- 

We will now explain how the multiplication by 
g,-(^-Eo)r jg carried out in continuous-time, that is with- 
out discretizing t. Consider an N-state system with 



Hamiltonian matrix elements: Hi. 



{i\H\j), where 




i,j = {l,2,..., TV}. For this A^-state system an instance 
of the (unnormalized) wave function is described by an 
A'^-dimensional vector having non-negative integer entries 



(2) 



This integer-valued vector represents M = Mi + M2 + 
. . . + Mpf replicas or copies of the system, where Mi of 
them are in state number 1, M2 are in state 2, etc.. For 
big systems the dimension of the vector is huge. With 
a finite number of replicas it will be very sparse and it 
is better to keep track of the state of each replica than 
writing down the vector explicitly. The requirement that 
the number of replicas in a given state is non-negative 
is rather restrictive, and is equivalent to the requirement 
that there is no sign problem. We will restrict oursclf to 
these cases. 

In order to build up the continuous-time formulation 
we will consider the evolution for an infinitesimal time 
step and then explain how to piece together (infinitely) 
many of these time steps in one shot. The action of the 
time evolution operator for an infinitesimal time step dr 
on an instance of the state is 

/ Ml 
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should note that a time varying Eb. causes the evolved 
wave function to deviate from the ground state wave 
bmctionfl^. This can be repaired by re-weighting the 
simulation as will be discussed in Sec. b 

We will now formulate a stochastic process that on 
average yields the evolution equation, Eq. ©. In the 
time interval dr a replica in state \i) can undergo one out 
of four different processes with associated probabilities: 

• "Transition", change state to j ^ i, probability 

• "Die", that is Mi ^ M, - 1, probability Poii)- 

• "Replicate", that is Mi Mi + 1, probability 

• "Stay", stay unchanged in state probability 

The "Die" and "Replicate" processes are known as 
branching processes. As these are all possibilities, prob- 
ability conservation implies 



(4) 



and must hold for all states i — 1, . . . ,N . 

The task of identifying the probabilities with matrix 
elements of the Hamiltonian is easy. Because the off- 
diagonal matrix element Hji is the only one responsible 
for transition between state i and j it is clear that 



(5) 



where j ^ i- In order to avoid the sign problem off- 
diagonal matrix elements are restricted to be negative. 

The increase in the number of replicas in state \i) from 
processes acting on replicas in state \i) is 



M' - M, 



Pr{i)-Pd{i)-Y.PtA^ 



M,. (6) 



This implies when comparing to the diagonal elements of 
Eq. (j3l and using Eq. © that 



PDi^)-PR{i)^ Hr. 




(7) 



(3) 

where the diagonal elements are Da = 1 + {Er — Hii)dT. 
Note that the ground state energy, £'0, is not known at 
the outset of the simulation. Therefore an estimator of 
the ground state energy known as the reference energy 
Er is introduced and used instead. During the course 
of the simulation this reference energy will be adjusted 
and can be used to extract the ground state energy. One 



The right hand side of the above takes either a positive 
or a negative value. In order to reduce the fluctuations 
/ in the replica numbers as much as possible Pr = is 
chosen whenever this value is positive and Pd = when 
it is negative. This choice implies that Po and Pr are 
of the order dr as also holds for Pt- The probability 
conservation equation Eq. then implies 



Psii) = l-\ \Hu-Ep 



dr, (8) 



FIG. 1: Selecting decay times according to the exponential 
distribution with decay constant A. Drawing a random num- 
ber r, the decay time Tdocay is selected as Tdccay = ln(r). 



which is of the order unity. 

We are now at the stage where we wish to patch to- 
gether many infinitesimal time steps. From the fact that 
Ps is of the order unity and all other processes are of 
the order dr it follows that for most (infinitesimal) time 
intervals nothing happens to a replica. The process can 
thus be modeled like the radioactive decay problem, al- 
though with several different decay channels: "Transi- 
tion", "Die", and "Replicate". The "Transition" channel 
is further divided into possibly as many as (N — 1) differ- 
ent channels corresponding to different non-zero values of 
Hji. This observation has been used previously to con- 
struct continuous-time algorithms for finite-temperature 
quantum Monte Carlo methods [20j . 

It follows that the imaginary-time evolution of one 
replica can be simulated by generating exponentially dis- 
tributed decay times with decay constant A = {Hu — 
Er + E,^. H,, I - H,., see Fig. □ 

Having obtained the decay time, the type of decay is 
determined stochastically proportional to the respective 
probabilities Pt,Pd and Pr, which are all of the same 
order dr. 

Practical implementation 

Although the method is formulated in continuous-time 
it is convenient to use discrete control times at regu- 
larly spaced time-intervals At as times where measure- 
ments are recorded and population control are being per- 
formed. 

In an actual simulation each replica contains informa- 
tion about the state of the system as well as a "clock" 
indicating the starting time for the next evolution. The 



FIG. 2; Illustration of a possible evolution for a 2-state 
system. The top line is the imaginary-time axis, on which 
control times are indicated by vertical bars. Initially there 
are three replicas, all in state 1 labeled by solid lines. As 
time progresses replica 1 first "Transitions" to state 2 (dotted 
line) , and then "Replicates" . Its replica also "Replicates" , but 
its copy "Dies" (the grounding symbol) almost immediately. 
Thus at the first control time the first replica has changed 
state and has divided into two. Nothing happens to the sec- 
ond replica up to the first control time, it just stays in the 
state 1. The third replica "Transitions" to state 2 just be- 
fore the first control time. Having propagated all three initial 
replicas and their descendants up to the first control time, the 
population control procedure is performed before propagation 
up to the next control time starts. 



replicas are ordered in a list. They all start in the same 
state and with their clock set to 0. Each replica in the 
list is subsequently evolved up to a control time Tc, or 
until the replica dies in which case it is removed from 
the list. An evolution of a replica begins by generating 
a decay-time Td according to the exponential distribu- 
tion. If Td > Tc the clock is set to Tc and evolution of 
the next replica in the list starts. If however t^ < Tc, the 
clock is set to Td and a random number is drawn to se- 
lect the decay type. If the decay type is "Transition" the 
state of the replica changes, if it is "Die" the replica is 
removed from the list, and finally the decay- type "Repli- 
cate" causes a copy of the replica with clock set to Td to 
be inserted at the end of the list. As long as the replica is 
not dead the evolution continues by picking a new decay 
time until Tc is reached. When the last replica in the list 
has evolved up to Tc, all replicas have the same clock- 
time, and measurements can be performed. The process 
is repeated by increasing and starting over from the 
beginning of the replica list. A graphical visualization of 
a possible evolution of replicas is shown in Fig. |21 

The control times are included in order to perform pop- 
ulation control to avoid an explosion/implosion in the 
number of replicas. Population control is achieved by 
changing the value of the reference energy Efj so as to 
maintain a roughly constant number of replicas. Ef; is 
only adjusted at the control times and is kept constant in 
between. Specifically, with E]^ denoting the reference en- 
ergy just after control time r^, which is the same value as 
just before the next control time r^"*"^, a possible choice 
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for population control is 



" " At \N^ 



(9) 



where Ni is the total number of replicas at and A^o 
is the total number of replicas at the beginning of the 
simulation. Eji denotes the average value of Ej^ for all 
control time points up to t^. Thus when the number 
of replicas decreases below Nq the reference energy is 
raised which will tend to increase the replication process, 
thereby increasing the replica population. 



Observables 

The measurement of observables in DMC requires 
some care 1,21:]. There are two issues which need to be 
addressed. First one must get rid of the dependence on 
the choice of initial and final wave function. Then one 
must be careful about getting rid of the bias introduced 
by the population control procedure. Lets first explain 
how to get rid of the dependence on the initial and fi- 
nal wave functions. This is known as the forward- or 
future- walking method j22|. 

When the initial wave function is a particular basis 
state \xi) the power method yields 



\xi) 



(10) 



where the overlap Cxi = (x/jV'o)- For simplicity in no- 
tation we have absorbed the reference energy into the 
Hamiltonian in this section. The projection can of course 
also be carried out for a wave function which is the super- 
position (with unit coefficients) of all orthonormal basis 
states \l) =Y.x 1^) 



(11) 



This can be converted into an evolution of the dual state 



(12) 



Note that the overlaps are all assumed to be real and 
positive (in accordance with the restriction of having no 
sign problem). 

In DMC the matrix evolution e^^'^\xi) is replaced by 
a stochastic process 



e-"^\xi) 



P{x,t; a;7,0)|a;) 



(13) 



where P(a;,r;x/,0) is the probability of finding the 
evolved state in the basis state \x) at time r provided 
it was in state \xi) at time zero. The stationary distribu- 
tion of this evolution is proportional to the ground state. 



It follows that the following can be used as an estimator 
of the ground state 



IV-o) 



N-r^OO 1 



CxrNr 

^ T X 

where the sum over different values of r contains Nr 
terms and the first r-value is taken after an initial equi- 
libration value Te- 

Using these relations one can obtain an expression for 
the ground-state matrix element of the observable de- 
scribed by the operator O 



r'.W^-tc 



aJ2iMe'"^'0\x)P{x, r; x/, 0) 



as well as for the norm of the ground state 



(15) 



•o) 



a^{l\e-"^'\x)P{x,T;xi,0). (16) 



To shorten notation we have collected the overlap coeffi- 
cients and the number of measurements into a single con- 
stant a = l/{CxjCiNr)- The overlap coefRcients cancel 
when considering the ratio 



(V'olV'o) 



(17) 



Diagonal observables 

Lets now specialize to the case where the observable 
O is a time-independent observable diagonal in the basis 
set Then the ground state matrix element is 

(V-olOl^o) 

= " E E(l|^') {^'\^^"''0\x)Pix, r; xj,0) 

T X.X' 

= aJ2J2^Mx')P{x', T + r'; x, t)OxP{x, r; x/, 0) 

T X.x' 

= aEE^(^''^ + ^''^'^)'^-^(^'^;^^'0) (18) 

T X.x' 

where we have used (1| = X]a;"(^"l' orthonormality of 
the set of basis states and the eigenvalue relation 012;) = 
\x)Ox- Similarly the norm is 

(V'oiV'o) = aEE^(^''^ + ^';^^'0) (19) 

T X' 

In both these expressions the limits r, t' — > oo are im- 
plied. This limit can of course not be achieved in practice 
so instead one evaluates the right hand sides with large 
values of r and r'. 
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Starting with typically thousands of replicas in the 
same initial state \xi) the probability P{x, r; xi,0) is pro- 
portional to the number of replicas in state \x) at time 
T. Thus the denominator is simply the total number of 
replicas at time r + r'. The numerator is a bit trickier 
in that one counts the replicas at time r + r' but with 
each replica weighted by the value of the observable at 
time T. To implement this in practice it is convenient to 
let each replica keep track of its history of observation 
values. This list must at least be of the length r'. Al- 
ternatively one could let each replica keep a history of 
its configurations, but this is quite memory-consuming, 
and it is better to concentrate on a few observables and 
their histories. The history list is such that whenever a 
walker replicates, the history list is also replicated such 
that the new walker inherits the same history as the orig- 
inal walker. 

The measurements must be performed at identical 
times for all replicas. One way of doing this is to measure 
the observable at the control times and store the result 
in the history list for each replica. One can also pick 
the measurement point to be at an arbitrary time in the 
time interval between replicas, instead of at precisely the 
control time itself. In fact, one can take the average of 
the observable over all time points in the interval, so as 
to get the maximum information available. This is quite 
easy to implement as the observable only changes value 
when a decay of the "Transition" type happens. As the 
times of decays are known it is quite easy to compute the 
time-weighted average of the observable for the interval 
between control times. For instance, assuming that there 
is just one decay at Td changing the observable from Oi 
to O2 , the average value accumulated at the control time 

is then (Oi(Trf - t^-^) + O^ir'^ " ^d))/(Ar). This av- 
erage value is then stored in the replica's history list of 
measurements . 



Off-diagonal observables 

It is slightly more complicated to measure off-diagonal 
static (time-independent) observables than diagonal ob- 
servables as there are no ways to insert off-diagonal 
operators without disturbing the configurations. The 
trick usually employed, which only works when the off- 
diagonal operator is a term in the Hamiltonian is best 
appreciated by visualizing the stochastic process as di- 
vided into small discrete time steps. That is 

P{x,t;xi,Q) = F(x,T;a;„_i,r„_i) . . . P(x2,T2\xi,ti)P{xi, 

(20) 

This discretization is never used in practice, it is used 
here merely for explaining the measurement procedure. 
In the limit of infinitesimal time intervals the processes 
involving changes of states are related to off-diagonal 
matrix elements of the Hamiltonian P{y, r + dr; x, t) = 



— (i/|_ff |a;)(iT, see Eq. The matrix element of the off- 
diagonal operator can then be written 



(V-olOlV-o) 



E 

x.x' ,x 



P{x',t' + t;x"t){x"\0\x) 
xP{x,T-dT;xi,0) (21) 



Note the "hole" introduced, there is nothing between r — 
dr and r. There is nothing wrong with having such a hole 
as the probabilities are invariant under time translations. 
Thus we can think of the left-most factor as being time- 
translated by an amount dr. Multiplying by 



1 



Po{x" , t;x,t — dr) 
-{x"\Ho\x)dT 



(22) 



one gets 



(Vo|0|V'o)= Pix',r' + T;x",T) 



y.Po{x" , t;x,t — dT)P{x, t — dr] xj, 0)- 



{x"\0\x) 



{x"\Ho\x)dT 
(23) 



where the subscript O indicates the transition described 
by the observable O. There is nothing special about 
the time r, other than it should be far from the start- 
ing and final time, thus we might take the average over 
many time intervals each of length dr within a range 
of nearby control times {tc~^,tI). Implementing this in 
practice amounts to counting the number of transitions 
corresponding to the off-diagonal operator in question 
within the time interval Ar and multiplying this num- 



ber by 



{^"\oW) 



This number is then stored as the 



-(x"\h\x')1:l7 

measurement value at time in the history list, and the 
variables are accumulated in the same way as for diago- 
nal observables using the forward- walking method where 
replicas are counted at t + t' weighted by the measure- 
ment value at r. 



Dynamic observables 

Dynamic observables can be recorded from the history 
list of measurement results for each replica. Dynamic cor- 
relation functions at times r = toAt significantly larger 
than the time interval between control times Ar can sim- 
ply be gotten by taking products of entries that differ 
rJzjycj^iO^lots in the history list. In order to measure the 
small time behavior using this approach Ar should be 
made smaller. Of course, as the method is formulated in 
continuous-time, there is nothing in principle restricting 
the measurement times to discrete time points mAr. At 
the expense of some extra book-keeping one can store 
and record observables at any time-separations. 
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Reweighting 

The use of population control where the reference en- 
ergy En is changed in order to keep the number of repli- 
cas approximately constant implies that the Monte Carlo 
procedure with population control is not sampling ex- 
actly the imaginary-time evolution of the wave function. 
Instead it is sampling an evolution with a time-dependent 
reference energy. Varying Ejj according to the recipe de- 
scribed in Sec. B the resulting wave function t/i will differ 
from the correct ground state wave function ^/^ by a prod- 



uct of time-dependent factors[2; 



^(r,) = e"SI=o(^«-^«)^"V(T-z) (24) 

The subtraction of the constant E'j^ is done in order to 
keep the exponential from overflowing. 

The extra multiplicative factor coming from the time- 
varying Eji can be gotten rid of by multiplying by the 
exponential factor above, giving for a diagonal observable 



(25) 



Of course one cannot keep track of an inflnite product of 
factors, and if one could, the fluctuations would be enor- 
mous, making it impossible to get accurate results [l^. 
However if the number of replicas isn't varying a lot, 
which can be achieved using importance sampling as de- 
scribed in the next section, the fluctuations need not be 

bigHi. 

In practice the observables are studied for different val- 
ues of reweighting times. Typically one sees a clear bias 
without reweighting which is possible to avoid by taking 
longer reweighting times, however too long reweighting 
times gives added noise. Fortunately the bias vanishes 
usually before the noise gets big, thus there is a region 
where one can get accurate measurements, see Fig. |3| 
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Importance Sampling 

It is known that importance sampling reduces statisti- 
cal errors in DMC[2^. Importance sampling is achieved 
by sampling the product of the wave function times a 
guiding function instead of the wave function alone. The 
guiding function is chosen to be as close to the exact 
ground state as possible. Lets now show that branch- 
ing can be avoided using importance sampling when the 
guiding function coincides with the exact ground state 
wave function. Consider the (infinitesimal time) evolu- 
tion equation 

(a|V(r -I- dr)> = (1 + {Er - H)dT) \a'){G'\i^{T)) 

(26) 

where \'iP{t)) labels the simulated instance of the wave 
function and |cr) denotes a basis state. Multiply by the 
time- independent guiding function {tjjg\o') and insert the 



FIG. 3: Columnar order parameter Xcoi on ^ 4x4 lattice vs. 
reweighting times for the quantum dimer model at V/ J = 0.1. 
M — 1000 replicas were used. The reweighting times are 
here defined using projections for r steps to project onto the 
ground state and forward-walking for r additional time steps. 
The line is the exact diagonalization result. 



unity factor ^^^j^,! 



(^,|a)(cT|^(T + dr)) = 

ii^.W) Y,{a\ (1 + (En - H)dr) |a')-|^(a'|^(r)). 



(27) 



As all eigenfunctions are real here, we do not distinguish 
between {ip\cr) and {cr\ip). The evolution matrix govern- 
ing the evolution of the product {i/jg\a){cr\il!{T)) is there- 
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fore 

{^g\a){a\{l + {ER-H)dT)\a')^^y (28) 

Now take the case where the guiding function ij^g coin- 
cides with the true ground state wave function. Then 
when we sum over columns of this matrix, that is 
and use the resolution of the identity, 1 = |cr)((T|, we 
get 

5](^g|a)(f7|(l + {Er - H)dT)\a')^-^^ 

^ {ijg\{En-Eo)dT\<j') 

(29) 

The last term vanishes provided Efj = Eq. The fact 
that each column sum to unity means that the matrix is 
Markovian. The evolution can then be simulated entirely 
without branching processes and is thus a classical Monte 
Carlo simulation in continuous-time. 

The diagonal elements of the evolution matrix are not 
affected by the guiding function, while the off-diagonal 
ones are. This changes the transitions between the states. 
Basically the transition to a state is enhanced if the guid- 
ing function has a large value for that state. If the state 
it was coming from also corresponds to a large guiding 
function value, tpg^ ensures that this is compensated for. 

QUANTUM DIMER MODEL ON THE SQUARE 
LATTICE 

We will now apply the continuous-time DMC method 
to perturbed RK-Hamiltonians taking the QDM on a 
square lattice as an example. The QDM was first pro- 
posed in the context of resonant valence bond (RVB) 
theories of high-Tc superconductivity fisj. In the RVB 
theory[2^ ji^l pairs of spins form singlets that are ap- 
proximated by dimers in the QDM. The Hamiltonian is 

H^~-jJ2 (|ll>(=l + H.C.) +VJ2 (ill) ("I + l=>(=l) 

(30) 

where the summations are taken over all elementary pla- 
quettes of the lattice. The basis states of the QDM con- 
sist of all possible dimer coverings of the lattice such that 
all sites are covered and no dimers overlap each other. 
The potential energy term V counts the number of pla- 
quettes possessing parallel dimers and the J-tcrm asso- 
ciates a kinetic energy to flips in the orientation of these 
parallel dimers. For this reason plaquettes possessing 
parallel dimers are named flipable plaquettes. 

Rokhsar and KivelsonQ realized that one could find 
the exact ground-state of this model a,t J = V. The 



ground-state at this RK-point is simply the equal super- 
position of all basis states 

10) = -^ El*) (31) 

i 

where \i) is a basis state describing a particular dimer 
covering of the lattice. N is the total number of basis 
states. Any ground-state expectation value of a time- 
independent operator diagonal in the basis states can 
then be evaluated as 

(0|^|0) = ^E^^ (32) 

i 

where the sum goes over all classical dimer covering 
states. The right hand side can be recognized as an ex- 
pectation value in classical statistical mechanics taken at 
infinite temperature, where the classical partition func- 
tion becomes Z — N. 

The special property of the QDM at = J which 
makes the simple wave function Eq. (|31|l the ground state 
is the fact that the Hamiltonian can be written as a sum 
weighted by positive coefficients of positive semi-definite 
Hermitean operators acting on pairs of states, where each 
operator is proportional to its own square. 

Eq. (|32|l describes a classical system at infinite tem- 
perature. The recipe for constructing finite-temperature 
RK-Hamiltonians was given in Ref. A quantum 

Hamiltonian of the form, which one can view as the gen- 
eral defining form of RK-Hamiltonians, 

+e''^'=^'-^^^/'\s'){s'\-\s){s'\-\s'){s\) 

(33) 

has the ground state 

l^o)-^E^"'''''^'l-) (34) 

with Eq = 0, where Z ~ e^^^" and is,s' res posi- 
tive. It follows that the ground state expectation value 
of static quantities can be gotten by evaluating a classical 
thermal expectation value at inverse temperature K. 

Quantum dynamics from classical Monte Carlo 

In addition to the "classical" ground-state properties 
of RK-Hamiltonians it was pointed out by Henley that 
the quantum dynamics of the RK-Hamiltonian can be 
gotten by performing a classical Monte Carlo simulation 
with the appropriate dynamics!^ |^. We want to show 
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that this relation between quantum dynamics and clas- 
sical Monte Carlo is general. This was also recently dis- 
cussed in Ref. 01 • Before doing so lets see what hap- 
pens to the continuous-time DMC method for the QDM 
exactly at the RK-point as well as for the generalized 
RK-Hamiltonians in Eq. 

A state with / flipable plaquettes in the QDM has 
potential energy fV. Having / flipable plaquettes implies 
that / other states are accessible by one dimer flip, thus 
there is a kinetic energy —fJ associated with this state. 
At the RK-point, V = J, the sum of potential and kinetic 
energy is thus zero. It follows from Eq. Q that Pjj — 
Pr = provided Eji is set equal to the true ground 
state energy Eq = 0. This holds for any basis state of 
the system. The fact that branching can be set to zero 
means that the DMC has been reduced to a classical 
Monte Carlo procedure. 

For the generalized RK-Hamiltonians described in 
Eq. H33|l , the branching processes are not automatically 
zero. However branching can be avoided if one simulates 
the product of the wave function with the ground state 
wave function Eq. H34|) . According to Eq. H28|l the Hamil- 
tonian then changes effectively, H — > ipgHtp~^, mean- 
ing that diagonal terms are unchanged while off-diagonal 
terms, Hs's e~^^^'''~^''-*/^i?s's, where we have used 
the ground state of the RK-Hamiltonian Eq. H34|l as guid- 
ing function. The amount of branching needed according 
to Eq. ^ becomes again 

- 0, (35) 

provided — Eq — 0. Thus the DMC algorithm from 
which quantum dynamics can be extracted reduces to 
a pure classical Monte Carlo simulation when used in 
combination with importance sampling using the exact 
ground state as guiding function. 

It is peculiar to note that in the quantum dimer model 
case the DMC method reduced to a classical Monte Carlo 
method without the explicit mentioning of a guiding 
function as was needed for the general RK-Hamiltonians. 
However, in fact the guiding function was used implicitly. 
This becomes clear when realizing that the ground-state 
Eq. (p?T|l is the equal amplitude superposition of all ba- 
sis states, thus '4ig{s)ijj~^{s') = 1 for all states s and s' . 
In other words the equal amplitude wave function is the 
default guiding function used in DMC when no other 
explicit guiding functions are specified. It corresponds 
simply to the multiplication by unity. 

Having explained how optimal importance sampling 
leads to a classical Monte Carlo procedure for RK- 
Hamiltonians it should be clear that the underlying rea- 
son follows from the known importance sampling argu- 



ment originally presented in Ref. [2J| and restated in 
Sec. Q; the operator governing the evolution of the wave 
function times the exact ground state wave function is a 
Markovian matrix. Thus quantum dynamics can be got- 
ten from classical Monte Carlo simulations for any Hamil- 
tonian that can be simulated with the DMC method pro- 
vided the exact ground-state is explicitly known. The 
advantage of dealing with RK-Hamiltonians is the fact 
that their ground states and energies are known explic- 
itly from the construction of the models. 

Measurement results 

Now turn to the physics of the QDM on the square 
lattice. For negative values of V it is favorable to have 
parallel arrangement of dimers around a plaquette. For 
positive V there is a potential energy cost to have such 
parallel arrangements, however the resonance term still 
favors parallel dimers. Thus it can be expected that the 
favorable dimer configurations maximize the number of 
parallel dimers also for positive V . In order to measure 
the dominance of parallel dimers we use the so called 
columnar order parameter which can be written as 

^coi = i^((E"H(^1(-l)'^^)' + (E"v(rl(-1)'^«)') 

(36) 

where nn (jt-v) is the number of horizontal (vertical) 
dimers belonging to the plaquette at r. r is an integer- 
valued coordinate labeling the center of each plaquette 
on a lattice of size N = L x L. The sums are to be taken 
over all plaquettes. 

The columnar order parameter is a diagonal observable 
which does not commute with the Hamiltonian. It is 
therefore necessary to use the forward walking technique, 
which require reweighting. 

Measuring the columnar order parameter for different 
system sizes at different values of V/J we obtain the re- 
sults shown in Fig. ^ |0| . It is seen that the columnar 
order parameter remains finite for V/ J < 1 and goes to 
zero atV/J=l. These results are consistent with earlier 
results obtained using exact diaeonalization fisl Eof . 

The columnar order parameter does not distinguish 
between columnar and plaquette phases, see Fig. |5l 

To distinguish between these it is helpful to consider 
the order parameter l3C| 

'^=^E("H(^1-"v(rO) (37) 

which is unity in the columnar phase and vanishes in the 
plaquette phase. A plot of (0^) as a function of V/J 
for different system sizes is found in Fig. is seen to 
decrease as V/J increases. This is consistent with the 
plaquette phase being more favorable as one approaches 
the RK-point. In order to detect a phase transition we 
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FIG. 4: Columnar order parameter Xcoi vs. V/J. The dif- 
ferent curves are for different linear system sizes L = 8, 16, 32 
and the open symbols are (quadratic) extrapolations to L = 
oo. 
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FIG. 5: Columnar phase (left) and plaquette phase (right). 
The presence of bonds indicates bigger probabilities for dimers 
than the absence of bonds. 



measure (0^) and {(j>*) and form the Binder cumulant 
5 = 1- (,/)4)/(3(0y) for different system sizes. The 
cumulants for different system sizes should cross at the 
phase transition provided the scahng regime has been 
reachedji^. Using exact diagonahzation the Binder cu- 
mulant was calculated for system sizes up to L = 8 us- 
ing and the phase transition was estimated to happen at 
roughly V/J = —0.2^^. While we reproduce the ex- 
act diagonahzation results for L = 4, 6,8 we find that 
for larger L the crossing points moves towards larger val- 
ues of V/J, see Fig. [3 It is rather difficult to pin down 
the exact location of the phase transition based on the 
crossing-points summarized in the inset of Fig. [7| How- 
ever it is clear that the estimate V/J = —0.2 gotten from 
exact diagonahzation studies is too low. Simulations on 
larger systems are needed in order to determine the ac- 
curate value. 

Off-diagonal observables can also be measured. Con- 
sider the operator F measuring the energy associated 
with plaquette ffips 

F = 5](|ll)(z|+H.c.) (38) 

where the sum is to be taken over all elementary pla- 




FIG. 6: The value of {(j)'^) vs. V/J for different system sizes. 




FIG. 7: Binder cumulant for the order parameter </!> vs. V/J 
for different system sizes. The upper right inset is a blow-up 
of the region containing the crossings of the curves for the 
largest systems sizes. The lower left inset shows the values of 
V/ J at intersections between L and L + 2 curves. 



quettes of the lattice. Figure |H1 shows the ground-state 
expectation value of F measured directly using Eq. H23|) 
(solid symbols) and indirectly measuring the potential 
energy which is a diagonal observable and subtracting 
that from the total energy (open symbols). It is evi- 
dent that for parameters where the guiding function is 
not optimal the direct measurement of the offdiagonal 
observable is more noisy than the indirect measurement. 

DMC has the advantage that imaginary-time correla- 
tion functions can be measured directly from the random 
walk. In Fig. 1^ we show the equal bond dimer-dimer dy- 
namic correlation function D — {Di{T)Di{Q)) — (Di)"^ 
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FIG. 8: Offdiagonal energy associated with plaquette flips vs. 
VjJ for a 16 X 16 lattice. The solid symbols denotes direct 
measurements of the offdiagonal terms and the open symbols 
are gotten from subtracting measurements of the potential 
energy {Eiocai) from the total energy (Etot)- 
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FIG. 9: semi-log plot of the equal bond dimer-dimer dynamic 
correlation function vs. imaginary time for different values 
of V/J. The curves are for a 16 x 16 lattice except for the 
two lowest curves which compares the Monte Carlo to exact 
diagonalization results on a 4 x 4 lattice. 



{Di is 1 if a dimer is present on bond i and otherwise). 
The lowest two curves show the agreement with exact 
diagonalizations on a 4 x 4-lattice (lines) for two values 
of V/J. Note the semi-log scale which "amplifies" the 
noise at small values of D. The upper two curves are 
for a 16 X 16-lattice. By fitting the long time behavior 
to an exponential we find that the finite size gaps for 
the 16 X 16-systems are 0.02 ± O.OIJ for V/J = 0.9 and 
0.022 ± O.OOIJ for V/J = 1.0. 



CONCLUSIONS 

We have explained in details a continuous-time formu- 
lation of the DMC method. The novelty of the method 
lies in its ability to output measurements of imaginary- 
time correlation functions without time discretization er- 
rors. 

As any DMC method it performs best when used in 
combination with a good guiding function resembling the 
ground state wave function. Good guiding functions ex- 
ist in the vicinity of RK-Hamiltonians, namely the ex- 
plicitly known ground state of the RK-Hamiltonian it- 
self. Thus the method performs well for perturbed RK- 
Hamiltonians. 

In the case when the guiding function is equal to the 
exact ground state the DMC method becomes equal to 
a classical Monte Carlo procedure. Thus in these cases 
imaginary-time quantum correlation functions can essen- 
tially be obtained by classical Monte Carlo. This was 
pointed out in Ref. 5] for RK-Hamiltonians, but does in 
fact hold generally whenever the ground-state is known 
explicitly. 

As examples of results that can be obtained using the 
continuous- time DMC we have calculated order parame- 
ters away from the RK-point in the QDM on the square 
lattice. 
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